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Abstract. Symmetry preserving difference schemes approximating second and third 
order ordinary differential equations are presented. They have the same three or 
four-dimensional symmetry groups as the original differential equations. The new 
difference schemes are tested as numerical methods. The obtained numerical solutions 
are shown to be much more accurate than those obtained by standard methods 
without an increase in cost. For an example involving a solution with a singularity 
in the integration region the symmetry preserving scheme, contrary to standard ones, 
provides solutions valid beyond the singular point. 
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1. Introduction 

The purpose of this article is to present some new difference schemes, having the same 
Lie point symmetry groups as the ordinary differential equations they approximate. We 
test these schemes as numerical methods and compare them with standard schemes. 

This is part of a general program, the aim of which is to turn Lie group theory 
into an efficient tool for solving difference equations. Continuous symmetries of discrete 
equations have been intensively studied during the last 20 years or so. 

For recent reviews containing extensive lists of references to the original papers, see 
[1, 2]. In this article we restrict ourselves to one specific aspect of this approach, the 
symmetry preserving discretization of ordinary difference equations and its applications 
in numerical analysis. 

To present the basic ideas, let us first consider an ordinary differential equation 
(ODE) 

E = yW-F(x,y,i/,...,y< n -V)=0. (1) 

Its Lie point symmetry group G consists of all local point transformations of the form 

x = A A (x, y), y = Sl x (x, y) (2) 

taking solutions y(x) into solutions y(x) of the same equation (A represents group 
parameters) . The Lie algebra L of the symmetry group G is realized by vector fields of 
the form 

X = £(x,y)dx + (f)(x,y)dy. (3) 

The algorithm for finding the symmetry algebra L and the symmetry group G for a 
given ODE (I) goes back to S. Lie and is given in many books on the subject [3]. It 
consists of solving the determining equations resulting from the infinitesimal invariance 
requirement 

W ^X(E)\ E=0 = 0, (4) 

where pr*™) X is the n-th order prolongation of the vector field X (acting on derivatives 
up to order n) [3]. 

Let us now consider an ordinary difference scheme (OAS), approximating the ODE 
(I). The scheme will consist of two equations relating the values of (x, y) in iV different 
points, with N >n + 1 

E a (n,x n+K ,. . . ,x n+L ,y n+K ,. . .y n+L ) = 0, a = l,2, L-K = N-1. (5) 

In the continuous limit one equation, say E 1 goes into the ODE (I), the other reduces 
to an identity (like = 0). The two equations (5) should be such that if AT — 1 values 
(xk,Uk) are given, we can calculate the iV-th one. This is assured e.g. by imposing a 
condition on the Jacobian: 
d(E 1 ,E 2 ) 



d(x n+L ,y n+L ) 



? 0. (6) 
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We wish to construct an OAS that not only approximates eq. (1), but has the 
same Lie point symmetry group G. This is achieved by constructing the scheme out 
of difference invariants of the group G, or out of invariant manifolds. These are found 
using the vector fields (3), corresponding to the invariance algebra of the ODE (1). The 
vector fields are the same as in the continuous case, however they must be prolonged to 
all points of the lattice, involved in the system (5). We have 

L 

pr X = ^2{£{x n+j , y n+j )dx n+j + <p{x n+j , y n+j )dy n+j }. (7) 

j=K 

The invariants satisfy 

pr XJ(x n+j , y n+j ) = 0, a = 1, . . . M, (8) 

where {Xl, . . . ,X M } is a basis of the algebra L. The invariant manifolds satisfy the 
same equation, but only on a subspace where the matrix of coefficients 

/ £l,n+K, ■ ■ ■ , £l,n+L, <f>l,n+K, ■ ■ ■ , 01,n+L \ 



\ £M,n+K, • • • , ^M,n+L, 4>M,n+K, • • • , 0M,n+L / 

is of lower rank. 

The fact that difference schemes can be invariant under continuous Lie point 
transformations that act on the equations and on lattices was pointed out by 
Dorodnitsyn [2, 4]. This approach has been used to classify and solve three-point 
difference schemes [5, 6]. It has been shown that symmetry preserving discretizations 
of first order ODEs are exact, i.e. the solutions of the ODEs and the invariant OAS 
coincide exactly [7]. The complementary problem, in which an OAS is given and we 
wish to find its Lie point symmetries was solved in [8]. 

In Section 2 we discretize a second order nonlinear ODE with a three-dimensional 
solvable symmetry algebra. Section 3 is devoted to a discretization of several third order 
nonlinear ODEs with three, or four-dimensional solvable symmetry algebras. Equations 
invariant under the simple group SL(2, R), or the reductive one GL(2, R) are discretized 
in Section 4. The obtained invariant difference schemes are tested in Section 5. They 
are shown to be considerably more accurate then the corresponding standard schemes. 
Moreover, in the study of a singular solution the invariant schemes turn out to have 
a qualitative advantage: they make it possible to integrate numerically beyond the 
singularity. 

2. Example 1 : A second order ODE invariant under a solvable Lie group 

Let us consider the second order ODE 

x 2 y" + Axy 1 + 2y = (2xy + xV) M/(fe_1) , k ^ 0, ±, 1, 2. (9) 

Its symmetry algebra has a basis given by 

9 2y d ^ r 1 d d ., . d , , 
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(for k — 0, | and 2 the symmetry algebra is larger and the equation is linear, or 
linearizable) . Equation (9) could be simplified by a transformation taking the algebra 
(10) into its standard form, but we are interested in discretizing it without prior 
simplifications. 

We note that the general solution of (9) is 

fc-i , 



y(x) 



k-1 



kx 2 



x 



xo) H 

x 



(11) 



where Xq and yo are integration constants. 

Now let us derive an OAS approximating the ODE (9), invariant under the Lie 
group generated by (10). We consider 3 points on a line x n _i, x n , x n+ i and the 
corresponding values y^ = y(x k ). The invariance condition 

prXF(x n - 1 ,x n ,x n+1 ,y n - 1 ,y n ,y n+1 ) =0 (12) 

with prX as in eq. (7) yields 3 elementary invariants: 



6 



2 2 



2 

X nVn 



X n-\Vn-\ 



(13) 



%n X n —\ {x n -\-\ X-fi) ' {x n X n —\) 

We put h n+ i = x n+ i — x n , h n = x n — x n -\ and expand y n ±\ = y(x n ±i) into Taylor series 
about x = x n . We obtain 



6 
l 

2 



(6 - TT-jfcZifc) = (/wi) 2_ *{ W + W + 2y) 



+ |(/^n+i - h n )(x 2 y"' + 6xy" + 6y') 



6 + 



(6) fc 



(fc-2)/(fc-l) 



= (^ t+1 ) 2 - fe (^V + 2a;|/) 



*(fc-2)/(fc-l) 



0(e 2 )} 
(14) 



x 



1 + ^(^1 



a; V' + 4zy' + 2y 
a; 2 ?/' + 2xy 



+ 0(£ s 



where we assume that h n+ i and h n are of order e. 
We see that the two equations 

- - I /--I — — •-. > ! = 



(6) 



1 (fc-2)/(fc-l) 



6 + (aF 16 



(15) 
(16) 



6 + 1 

with = const, provide an invariant OAS approximating eq. (9). In general this is a 
first order approximation (of order e). If we choose K = 1 in (16), the first order terms 
drop out and we obtain a second order scheme (and a uniform lattice). 



3. Examples of third order ODEs invariant under solvable Lie groups 

3.1. General comments 

A third order ODE can have a Lie point symmetry group of dimension dim£ = N, 
< N < 7. The maximal dimension N = 7 occurs only for linear equations that can be 
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transformed into y'" = by a point transformation [3, 10]. We shall consider examples 
of equations with N = 3 and N = 4. 

In order to approximate a third order ODE we must consider at least 4 points in a 
stencil. We denote the points 

(x n+k ,y n+k ), -1 < k < 2, (17) 

and put 

h n +2 = x n+2 ~~ x n+l) = x n+l ~~ x ni h n = X n — X n -\. (18) 

In the continuous limit we put 

h n+j = ctjE, j = 0,1, 2 (19) 
where otj are constants of the order of 1 (not necessarily all equal). 

3.2. Example 2 : An ODE invariant under the similitude group of a Euclidean plane 

Let us consider the four- dimensional Lie algebra 



Ox Oy Ox Oy Ox Oy 

generating translations, rotations and dilations in the (x, y) plane, respectively. The 
Euclidean algebra {Xi, X 2 , X 3 } allows two independent differential invariants in the 
space {x, y, y', y", y'"} 

T = y" r = (l + y' 2 )y'" - Zy'y" 2 ( . 

1 (l + y' 2 ) 3 / 2 ' 2 {l + y' 2 ) 3 { ' 

and the invariant ODE is 

h = F(h), (22) 

where F(z) is an arbitrary function. 

Invariance under dilations corresponding to X 4 implies F(z) = Kz 2 and the 
invariant ODE is 

(1 + y ' 2 )y"> - 3y'y" 2 = Ky" 2 (23) 
where K is a constant. The general solution of eq. (23) can be given in implicit form as 

u(t)dt + C 3 , x = C 1 J^^—^ds + C 2 , (24) 

where Ci, C 2 , and C3 are constants. 

The Euclidean Lie group corresponding to {Xi, X 2 , X 3 } allows 5 functionally 
independent difference invariants in the space with local coordinates (17). We choose 
the following basis for the invariants: 

, / Vn+2 ~ Vn+1 
?1 — lln+2 J- 



h n +2 

2-1 1/2 

C2 = h n+ \ 



-y _|_ / Vn+1 — Vn 



h 



n+1 
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2-1 1/2 



6 = K 



Vn - Vn-1 



6 

(25) 



£4 = (y n +2 - y n +i)h n+1 - (y n+1 - y n )h n +2, 
6 = (y n +i - y n )K - (y n - y n -i)h n+1 . 

From these we can form invariants that approximate the differential invariants l\ and 
h of eq. (21). To see this we expand y n+2 = y{x n + h n+l + h n+2 ), y n +i = y(x n + h n+1 ) 
and y n -i = y(x n — h n ) into Taylor series about x n and obtain 

6/6 6 



J9 



1 + y' 2 )y>" - 3y'y 



1 • '2 



6 + 6 + 6 + 66(6 + 6) 

1 

" (1 + y' 2 ) 3 

^-n+2 + 2/l n+ i — Al n 



(26) 



+ 



1 + y' 2 )y™ - I0y'y"y"' + 15y 2 y 



2 "3 



3 y'" 3 2h 2 n+2 + 7h n+2 h n+1 +4h 2 n+1 + h n+1 h n -2h 2 r 



8 1 + y' 2 



{h n+ 2 + h n+ \ + h 

n ) 



Jl = 



2«6 



+ 



2/^6 



66(6 + 6) (66X6 + 6) 

(i + y^ ^ + ^-^[(1 + y' 2 )y'" - zy'y"^ 



(27) 



3(1 + y' 2 ) 

x [a(h n+2 + 2h n+1 ) + /3(h n+1 - h n ] j , 

a + /3 = 1. 

An invariant OAS approximating eq. (22) is given by 

J2 = F(J ± ) (28) 

on the lattice: 

Mi + ^6 + C6 = 0, (29) 

where A, B and C are constants. 

In particular the ODE (23) invariant under the similitude group Sim(2) is 
approximated by 

J 2 = KJ 2 (30) 

on the lattice (29) which is also invariant under Sim(2). 

Other invariant lattices can be formed out of the invariants (25), for instance in 
Subection 5.3 we choose 

6_6 
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Generally speaking (29) and (30) (or (31)) provide a first order approximation (of 
order e if h n+2 , h n+1 and h n are of order e). For a special value of K we can cancel first 
order terms in e and obtain a second order approximation, namely 



K 



a 



2> 



C=-A, B = 2A. 



3.3. Example 3 : Equations invariant under a Euclidean Lie group 

Let us consider a different realization of the Euclidean and similitude Lie algebras, 
namely 



d 



x- 



d 



d 



x d 



y- 



d 



(32) 



d 

— , X 2 = x— , X 3 = (1 + x z )— + xy- ^ „ 
ay ay ox ay ay 

This algebra is isomorphic to (20) but cannot be transformed into it by a transformation 

of variables. The Euclidean Lie group corresponding to {Xx^Xi^X^} allows two 

independent differential invariants of order 3 or less. We choose them to be 

h = (i + x Y 2 y", (33) 
I 2 = [(l + x ^y'" + 3xy"}(l + x 2 ) 3 / 2 - 
The invariant third order ODE is 

h = F(h) (34) 

where F(z) is an arbitrary function. If we also require invariance under the dilations 
generated by X 4 , we obtain F(z) = Az and the equation is a linear one. 

The five functionally independent difference invariants in the space (17) allowed by 
the Euclidean group generated by {X-y, X 2 , X 3 } are 

Vn+l ~ Vn Dn ~ Dn-l 



6 
6 



(1+^) 1/2 

(l+x 



Xn+1 



Xr 



Xr 



Xn—1 



X r 



2 \l/2 
n+l) 

%n—l 



y n +2 - y n +i y n +i - y n 



x n+2 ~ x n+l 



X. 



n+l 



X r 



(35) 



x 



n+l 



c x n+2 ~ x n+l 
Q — T— , ?4 — -T— , ?5 — T— • 

1 -+- X n X n _\ 1 -+■ X n X n+ i 1 -+■ X n+ iX n+ 2 

Expanding into Taylor series about the point x = x n we find 



26 



6+6 

26 

6+6 



[i+x 2 f 2 

{l+x 2 f 2 



y"+ 



h n +i~h r . 
3 

h n +2~ 



/// Od, n n 



l+x: 



-y 



2h, 



y'"+ 



i+xr j 



+ 0(e 2 ) 
"] +0(e 2 ) 



(36) 



We have assumed that h n , h n+ i and h n+2 are all of order e (but not necessarily equal). 
From eq. (36) we obtain 



Ji = 
J 2 = 



2a6 + 2/36 



6 + 6 
6 



6 + 6 



ll + x 2 ) 3 ' 2 y" + 0(e), a + P = l, 



(37) 



6 + 6 + 6 V6 + 6 



6 



6 



6 + 6, 

= (l + a; 2 ) 3 / 2 [(l + a;V 



3xy"} + 0{e). 
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Thus an invariant OAS approximating eq. (34) is given by 

J 2 = F(J l ), 0^3 + ^4 + ^5 = (38) 

where a, b and c are constants. In general this will be a first order scheme. For certain 
functions F the scheme can be improved to a second order one by an appropriate choice 
of the constants a, (3, a, b and c. We shall not go into that here. 



4. Third order equations invariant under SL(2,R) 

Four inequivalent realizations of sl(2,R) as subalgebras of diff (2, K.) exist [5, 11]. We 
shall consider two of them here. 



4-1. Example 4 ■ First sl(2,R) algebra 
The first si (2, R) algebra Si is given by 



X i = 7T' X 2 = + y-^-, X 3 = 2xy-^- + y 2 -^-. (39) 
ay ox ay ox oy 

It allows one second order and one third order differential invariant: 

T 2xy" + y' T x 2 (y'y"' -3y" 2 ) 

y 3 y 5 

The most general invariant third order ODE is hence 

h = F(I 1 ), (41) 

where F(z) is an arbitrary function. The algebra (39) can be extended to gl(2,R) by 
adding the operator 

X 4 = *|. (42) 

Requiring invariance under the corresponding GL(2, R) group restricts F(z) to F(z) = 
Az 3/2 and we obtain the ODE 

x 2 (y'y"> - 3y" 2 ) = Ay'^xy" + y'f 2 . (43) 

Five independent SL(2,R) difference invariants are 

-.{Vn+2 ~ Vn+l), 6 = ; 1 (Vn+l ~ Vn), 



\J %n+l%n+2 \J %n%n+l 

1 1 
6 = ■ (Vn - Vn-l), 6 = (y n +2 - Vn), (44) 

6 = (Vn+l ~ Vn-l)- 

V •En+l%n—l 



From these we form 



T = (6 - 6 - 6X6 + 6)6 - (6 - 6 - 6)6(6 + 6) , 4 . 

666(6 + 6)(6 + 6)(6 + 6 + 6) 1 J 

x 2 (y>y'"-3y" 2 ) 
= r, V £<P, 

y 5 
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Ji = 8 



£4 — £l — £2 /-. _ s £5 — £2 — £3 



(46) 



2xy" + y' 2 x(y'y"' — 3y 



'2 



, +-[a(^ ++ + 2/i + ) + (l-a)(/i+-/i)]: 
2/ 1/ 



where and ^ are some functions of re, y', y", y'" and and < a < 1 is a constant. 
An invariant scheme approximating eq. (41) is given by 

J2 = F{J 1 ), (47) 
Mi + ^6 + C£ 3 + DU + ££ 5 = 0, (48) 
where A, . . . ,E are constants. To lowest orders (48) yields: 

v' 

{Ah n+2 + Bh n+1 + Ch n + D(h n+1 + h n+2 ) + E{h n+1 + h n )}?- (49) 

x 

+ {A(h n+2 + 2h n+1 )h n+1 + Bh 2 n - Ch 2 n _ x + D(h n+1 + /i„ +2 ) 2 

In general, the scheme is a first order one, i.e. all hj go to zero like hj = OjS, then the 
error in (47), (48) goes to zero like e 1 . For specific functions F(z) the accuracy can be 
improved by an appropriate choice of the constants a and A, . . . , E. 

4-2. Example 5 : Second sl(2, R) algebra 

The second sl(2,R) algebra S 2 has a basis given by 

= X 2 = |/|-, X 3 = y 2 ^-. (50) 

ay ay ay 

It can be embedded into the algebra sl(2, R) © sl(2, R) by adding 

X 4 = — , X 5 = a;— , X 6 = x 2 — • (51) 
or ax ax 

The Lie group SL(2, R) generated by S 2 has two differential invariants in the considered 

space, namely 

h = K (y T - h" 2 ), h = x. (52) 



y >2 y * 2 

The invariant ODE is 

^(yT-\y' 2 )=F{x). (53) 

Requiring invariance under the GL(2, R) group that includes X 4 in its Lie algebra 
reduces (53) to 

^(yy"-l y ' 2 )= K (54) 

where K is a constant. 
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A larger invariance group is obtained only for K = 0. In this case the equation is 
invariant under SL(2,R) <g> SL(2,R), generated by (50) and (51). 
The difference invariants corresponding to the algebra (50) are 

(y n +2 - y n )(y n +i - y n -i) h , , 

IX—-, —, r, X, lL n +2i "■n+li ,l n- 1<-> J ) 

(yn+2 - Vn+lAVn ~ Vn-l) 



We have 



6h n+2 h n 

h n +i(h n+ i + h n+2 ){h n + h n+ i)(h n+ 2 + h n+ i + h n ) 



x 



(h n+ 2 + h n+ i)(h n+ i + h n ) _ ^ 
h n h n +\ 



'0 

y 



I III 3 "o 

y y - 



+ 0(e). 



An invariant OAS approximating eq. (53) is 

J\ — F(x n , h n , h n+ i, h n+ 2), <f>(x n , h n , h n +i, h n+2 ) = (57) 

with 

F(x n ,0,0,0) = F(x), (58) 
<j>(x n , 0,0,0) =0. (59) 

If we require invariance under the group corresponding to {X±, X 2 , X 3 , X 4 } we must 
take F(x) = K and the lattice will depend only on h n+2 , h n+ i, and h n . For instance we 
can take the lattice to be given by 

ah n+2 + f3h n+1 + jh n = (60) 

and the constants a, (3, 7 can be chosen to improve the approximation. 

An OAS invariant under SL(2, R) <g> SL(2, R) that approximates eq. (54) for K = 

is 

(x n+2 - x n ) (x n+1 - x n -i) {y n+2 - y n ) (y n +i - y n -i) 



(x n+2 - x n+l )(x n - x„_i) (y n+2 - y n +i){y n - y n -i) 

(%n+2 ~ %n) (^n+l — x n-l) 



= 0, (61) 



„ +1 ) = ^" m 

For K = 4 this scheme is an exact one. Indeed, the equation 

y'y'" - \y" 2 = o (63) 

has two families of solutions 

y = + c and y = ax + (3 (64) 

ax + b 

where a, b, c, a and (5 are integration constants eq. (62) with K Q = has two families 
of solutions 

x n = + c, x n — an + (3. (65) 

an + b 

On the lattice (65) the functions (64) solve (61) exactly. 
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In this example the underlying Lie group SL(2, R) plays a specially prominent role. 
The group is the group of projective transformations of the real line (the variable y). Its 
invariant l\ is the Schwarzian derivative of the variable y. Projective transformations 
can be used to transform any three points on the projective line into any other three 
chosen points. Given four points, e.g. y n -i, y n , Un+i, Un+2 we can form precisely one 
projective invariant out of them, namely the anharmonic ratio R of (55). The variables 
x n , h n+2 , h n+ i, h n in (55) are also invariants since the considered SL(2,R) group acts on 
the y space only. We can call (53) a Schwarzian ODE. Then (57) is a Schwarzian OAS. 
Schwarzian derivatives play a prominent role in the theory of integrable systems [16] 
and of dynamical systems [17, 18]. 

5. Numerical results 

5.1. General procedure for testing the numerical schemes 

This section reports on the numerical experiments performed using the schemes 
described in the previous sections. The schemes are used to compute the solution 
for initial value problems on a given interval. Before describing the results for each of 
the four classes of symmetries analyzed in this paper, we first describe some general 
procedures to implement and test the various methods. 

5.1.1. Reference solution For test-problems for which an analytical solution is not 
available, a very accurate and reliable reference solution is computed numerically and 
used to assess the performance of the point symmetry preserving scheme. This is 
done using Matlab's standard adaptive Runge-Kutta scheme ODE45, with a very strict 
tolerance on the error set at tol = 10~ 9 . The first step is to convert the n-th order 
Equation (1) for y(x) into a system of n first order ODEs for u\(x) = y(x), u 2 (x) = y'(x), 
. . . , u n (x) = y^'^^x). Then Equation (1) becomes the system 

u[ = u 2 , u' 2 = u 3 , u' n = F(x, ui, u 2 , ■ ■ ■ , u n ). (66) 

Given initial conditions Ui(x ), u 2 (x ), . . . ,u n (x ), one then proceeds to compute the 
solution on the interval [xq, xf], where the scheme adaptively selects the local integration 
step so that its local error estimates satisfies the imposed tolerance. Those very high 
order, very accurate (and very costly numerically) solutions are used to generate start- 
up values as well as error estimations for the point symmetry preserving schemes as 
described next. 

5.1.2. Start-up values The symmetry preserving schemes require a number of start-up 
values (y = y(x ), y\ = y(xi) for the second order case; also y 2 = y(x 2 ) for the third 
order cases). For given initial values y(x ), y'(x ), (and y"(x ) for the third order case), 
the start-up value yo = y(xo) is directly available, while the values for y\ and y 2 are 
obtained as the the numerical reference solution (obtained as described above) at the 
nodes x\ and x 2 . 
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5.1.3. Error analysis Given the discrete mesh x n , n — 0, 1, 2, . . . , TV and corresponding 
solution y n generated by the point-symmetry preserving scheme, the corresponding 
errors are obtained by comparing y n with y re f{x n ). Although the user has no direct 
input on the actual mesh used by the Matlab's solver, it is possible for the user to 
request specific output points for the discrete solutions, so that given x n , one can obtain 
a very reliable numerical approximation y Te {(x n ), accurate with the prescribed tolerance. 

5.1.4. Equivalent standard schemes To better assess the new schemes proposed here, 
their performance for various test-cases is compared with that of the standard finite 
difference schemes that uses the same number of grid points x SjTl , n = 0, 1, 2, . . . , N. 
Although the point-symmetry preserving scheme finite mesh is typically non-uniform, 
for simplicity, the standard mesh is assumed to be, so that x s>n = Xq + nh with 
h = (xf — xq)/N. The discrete standard scheme is obtained using the following standard 
procedure (given here for the third order case, easily adapted for the second order case). 
Given the four points (z s ,„-i, y s , n -i), (x s ,n,y s ,n), (x s ,n+i,y s ,n+i), (x s ,n+2,y s ,n+2)- 

(i) obtain the interpolating polynomial Ps(x) through the four given points 

(ii) evaluate analytically i^(x a>n+ i/ 2 ), P3 (z s , n+ i/ 2 ), P^{x s ^ n+1/2 ), which gives: 

^3(2^+1/2) = ^-(27(y s , n+ i - y s>n ) - (y s , n +2 - y s ,n-i)) (67) 

Pz{x S ,n+l/2) = 2^(y S ,n+2 ~ (Vs.n+l + Vs,n) + Vs.n-l) (68) 

P^{x Stn+ i /2 ) = -^{y s , n +2 - 3y s , n +i + 3y s , n - y s ,n-i) (69) 

(iii) substitute those expressions in the equation being discretized, evaluated at x — 

X s ,n+l/2- 

5.2. Numerical experiments for Example 1 (second order ODE (9) ) 

Selecting k = 3 in Equation (9) gives: 

x 2 y" + 4xy' + 2y = (2xy + x 2 y') 1/2 (70) 

to be solved for x in the interval [1, 3], with the initial conditions chosen as y(l) = 13/12, 
y'(l) = —1. The problem has an exact solution: y r ef(^) — x/12 + 1/x 2 . 

The symmetry preserving scheme (15), (16) is used with the special choice K — 1, 
so that the mesh is uniform and the discrete scheme is given by: 

x 2 n+1 y n +i ~ 2x 2 n y n + x 2 n _ x y n ^ = (|) 1/2 /i 3/2 (a; T 2 t+1 i/ n+1 - x* n _ x y n -i) 112 (71) 

with x n = Xq + nh. 

The start-up values are given by y = y(x = 1) = 13/12 and y\ = y(x = 1 + h) = 
(l + /i)/12 + l/(l + /i) 2 . The corresponding standard discrete scheme is given by: 

(y s ,n+i - 2y s ,n + y s , n -i)x\ n + 2x S)n h(y s 

,n+l ys,n— 1 

) + 2h 2 y s>n (72) 

/ \ 1/2 

_ 7 2 / r> _|_ 2 y^n+1 ys,n—l \ 

— n \zx S!n y Stn -+- x s n — i 
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Table 1. Discretization errors, Example 1 



Scheme 


h = 0.1 


h = 0.01 


h = 0.001 


Sym.pres. 
Standard 


6.04 10" 4 
4.72 10~ 3 


7.26 10- 6 
7.54 10~ 5 


7.39 10- 8 
7.86 10~ 7 



Note that both the symmetry preserving scheme as well as the standard scheme lead 
to a nonlinear problem to compute y n+ i, given y n and y n -i- Those algebraic nonlinear 
problems are solved using a standard fixed point iteration until convergence. 

Using the exact solution as reference, errors are computed for the numerical 
solutions using each of the two schemes, with mesh sizes h = 0.1, 0.01, 0.001. Those 
errors are reported in Table 1. 

One observes from Table 1 that both schemes are second order accurate, as the error 
is roughly divided by a factor 100 whenever the mesh size is divided by 10. Also, the 
errors from the symmetry preserving schemes are smaller by a factor of 10 compared 
to the errors obtained with the standard scheme with the same mesh size. This is 
achieved without any additional computational cost, both schemes having the same 
computational complexity. Figure 1 shows the error as a function of x for both schemes 
for mesh size h = 0.1. The gain from using the symmetry preserving scheme is obvious. 



5. 3. Numerical experiments for Example 2 ( third order ODE (23) ) 

The test-case consists in solving Equation (23) for K — 1, with x in the interval [0, 10] 
and with initial values y(0) = 0, y'(0) = —10, y"(0) = 1. The lattice equation is chosen 
in the form (31), i.e. 

-r = T = 1 - (73) 

Start-up values for xo — 0, x\ — ho, X2 = 2/i are generated for a given ho using 
the Matlab solver, see Subsection 5.1. The constant 7 in (73) is then computed using 
the start-up points (x ,y ), (x 1 ,y 1 ), and (x 2 ,y 2 )- 

Given the three points (x n -i, y n -i), (x n , y n ), (x n +i,y n +i), the new point (x n+2 , y n +2) 
is obtained as the solution of the nonlinear system consisting of Equations (30) (with 
K — 1) and (73). In the present set of experiments, the values a — (5 — \ were 
selected. The resulting problem for (x n+2 ,y n +2) is nonlinear, in particular the mesh x n 
is non-uniform and completely coupled with the solution y n . 

The standard scheme is obtained by substituting the expressions in (67) (68) (69) in 
(23). It also leads to a nonlinear problem for y n+ 2, but for that scheme, the mesh x n is 
assumed to be uniform and certainly completely decoupled from the solution y n . 

Table 2 reports the numerical errors corresponding to various values for ho to start 
up the symmetry preserving schemes: ho = 1, 0.1, 0.01, which lead to respectively 14, 
130 and 1297 mesh nodes. The solutions with the standard schemes were computed on 
uniform meshes with the same number of nodes. 
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Table 2. Discretization errors, Example 2 



Scheme 


h = l 
(N = 14) 


h = 0.1 
(N = 130) 


h = 0.01 
(N = 1297) 


Sym.pres. 
Standard 


2.14 1(T 5 
4.20 10~ 2 


2.98 10~ 7 
5.83 10" 4 


6.45 10~ 9 
6.01 10" 6 



Both schemes appear to be effectively second order accurate, with the error in the 
symmetry preserving scheme smaller by a factor 1000. The discretization errors with 
both schemes are shown in Figure 2. 



5.4- Numerical experiments for Example 3 (third order ODE (34) ) 

The test-case consists of solving Equation (34) for the special choice F(Ii) = If, which 
leads to the equation: 

(1 + x 2 )y"' + 3xy" = y"\\ + x 2 f /2 (74) 

The solution is sought for x in the interval [0, L], with L to be selected below. The 
start-up values (x = 0, y = y(x )), (x 1 = h ,yi = y(zi)), (x 2 = 2h ,y 2 = y(x 2 )) are 
obtained as before using an over- resolved numerical integrator. The procedure to generate 
the mesh x n , n = 0, 1, 2, . . . , N and the corresponding discrete solution y n is as follows: 

• Step 1. Using the invariant Equation (38), one generates the complete mesh (for 
this particular case, it is independent of y n ). The constants in Equation (38) are 
taken as a = 1, b = —7, c = leading to 

T = T =1- (75) 

The strategy to select 7 and compute the corresponding mesh is the same as the 
one used for Example 2, see discussion above. 

• Step 2. Given the mesh x n , solve the invariant Equation (38) for y n+2 given 

(Xn-l, Un-l)-, {XniUn)-, (^ra+1 , Vn+l) , and X n+2 . 

Noting that I 2 = (1 +x 2 Y 1 ^dIi/dx, the equation being solved can be rewritten as 
(1 + x 2 )dli/dx = If. The solution for h(x) is therefore given by: 

— = arctan(a;) (76) 

h h,o 

where I±fi = h(xo = 0). This shows that y"(x) will blow-up if x — tan(l//i i0 ). We 
assess the performance of the scheme for two cases, one with blow-up and one without. 

5.4-1- Blow-up case The integration is performed for x in the interval [0, 11.2] with 
blow-up set up to occur at Xb = 11.25. This is achieved by imposing y"(0) = 
1/ arctan^b). Three values for the initial ho are selected to be ho = 0.1, 0.01, 0.001, 
which lead to meshes with respectively iV = 18, 151, 1484 nodes. Table 3 reports the 



Difference schemes with point symmetries and their numerical tests 15 
Table 3. Discretization errors, Example 3, case with blow-up. 



Scheme 


h = 0.1 


h = 0.01 


h = 0.001 


Sym.pres. 
Standard 


8.82 10~ 2 
7.03 10- 1 


9.88 10- 3 

1.10 lo- 1 


3.93 10~ 4 
1.68 10~ 3 



Table 4. Discretization errors, Example 3, case without blow-up. 



Scheme 


ho = 0.1 


h = 0.01 


h = 0.001 


Sym.pres. 
Standard 


1.53 10~ 3 
1.44 1Q- 1 


1.62 10~ 5 
9.65 10~ 3 


2.63 lO" 6 
1.18 10~ 4 



errors at Xp = 11.2. Both the symmetry-preserving and the standard schemes appear 
to be of order 1, with the errors from the symmetry-preserving schemes significantly 
smaller. 

Figure 3 shows the behaviour of the discretization errors for both schemes, for the 
case h = 0.01. 



5.4-2. No blow-up case This time, we select y"(0) = —1/ arctan(a;b), so that blow-up 
will not occur for x > 0. The numerical experiments are repeated with this new initial 
value. Table 4 presents the errors at xf for various values of ho, the conclusions are the 
same as for the blow-up case: both schemes appear to be first order accurate, with the 
symmetry-preserving scheme much more accurate. 
Figure 4 illustrates this behaviour for h = 0.01. 

5.5. Numerical experiments for Example 4 (third order ODE (43) ) 
The test-case consists of solving Equation (43) with A = — 1, i.e. the difference equation 
h = -JT (77) 
on the lattice given by 

j- = 7- (78) 

with J 2 , J\ as in (45) and (46) and as in (44). The solution is sought for x in the 
interval [1, 16] with initial conditions y(l) = 0, y'(l) = 0.1, y"(l) = 0.1. The start-up 
values (x = l,y = y(x )), (xi = 1 + h ,yi = y(xi)), (x 2 = 1 + 2h ,y 2 = y{x 2 )) 
are computed as in the other cases. Given the three points (x n _i, y n -i), (x n ,y n ), 
(x n+ i,y n+ i), the next point (x n+ 2,y n +2) is obtained as the solution of the nonlinear 
system corresponding to the two symmetry preserving discrete Equations (77) and 
(78). The constant 7 in (78) is computed based on the three start-up values. Table 5 
contains the errors with the symmetry preserving scheme and the standard scheme for 
this example, corresponding to various initial mesh sizes ho = 0.2, 0.01, 0.005. 
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Table 5. Discretization errors, Example 4 



Scheme 


h = 0.02 
(N = 149) 


h = 0.01 

(N = 288) 


h = 0.005 
(N = 567) 


Sym.pres. 
Standard 


2.67 10~ 4 
1.47 10- 3 


6.62 10~ 5 
5.59 10" 4 


1.65 10~ 5 
1.78 10~ 4 



According to the results in Table 5, both schemes are second order accurate, with 
a much smaller error for the symmetry preserving scheme. Figure 5 represents the 
discretization error behaviour for h = 0.04. 



5.6. Numerical experiments for Example 5 (third order ODE (53) ) 

Numerical experiments are conducted with Equation (53) for the case F(x) = sin(x): 

^Y'-\y 2 ^sin{x). (79) 

The solution is sought for x in the interval [0, 2] (also [0, 6]) with initial conditions 
2/(1) = 0, y'(l) = —10, y"(l) — 1. A uniform mesh is used here, it is compatible with 
the difference invariants in (57). With h n = h n+ i = h n+2 = h corresponding to the 
uniform mesh, the other invariant difference equation in (57) becomes: 

R = A(l-^-F(x n ,h)) (80) 

with R defined in (55) as R = (y n+2 - y n ){Vn+i - Vn-i)/ ((Vn+2 - yn+i)(y n - y n -i)) and 
where we select F(x n ,h) = F(x n + h/2) to achieve second order accuracy. This leads 
to the following explicit expression for y n+ 2- 

(y n +i - Vn-\)Vn ~ K(y n - y n _i)y n+ i 

Vn+2 = s TT-, s (.81 J 

{Vn+l - Vn-i) - K(y n - y n _!) 

where K = 4(1 — ^-F(x n + h/2)). This explicit expression for y n +2 is remarkably simple. 

On the other hand, the standard scheme for the same problem is nonlinear. 
Substituting the finite difference approximations for y',y".y"' in (67), (68), (69) in the 
ODE (79) leads to a nonlinear equation for y n+ 2 to be solved iteratively. 

First, we compare the discretization errors using the invariant scheme and the 
standard scheme on the interval [0, 2] on which the solution is smooth. Table 6 reports 
those errors in terms of the mesh size h. Both schemes display a second order convergence 
rate. The standard scheme has errors which are smaller by a factor of 6, but in terms of 
computational efforts, the invariant scheme is much more efficient, as it gives an explicit 
formula for y n +2 unlike the standard scheme that requires a nonlinear iterative solver at 
each step. However, if the integration interval is [0, 6], remarkably different conclusions 
are obtained. The solution develops a singularity around x — 3. At that point, both 
the standard scheme and the adaptive Runge-Kutta solver from Matlab fail to converge. 
On the other hand, the invariant scheme integrates right through the singularity. The 
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solution obtained with the three schemes (reference, standard, invariant) is displayed 
in Figure 6 for the coarse resolution h — 0.1. In Figure 7, the solution is shown with 
the invariant scheme for three resolutions h = 0.1,0.01,0.001. To better observe the 
behavior of the solution near the singularity, the plot uses a log scale, and the absolute 
value of the solution is shown. Excellent numerical convergence is observed, with the 
solutions corresponding to the three resolutions matching very closely each other (of 
course, the singularity is better captured by the finest mesh). 

The most striking feature shown on Figure 6 and 7 is that the symmetry preserving 
difference scheme provides a numerical solution u(x) for the entire region < x < 6 , 
x xq, even though the solution has a pole at xq close to 3. A similar phenomenon 
was observed in a previous study of a specific type of first order systems of ODEs, 
namely matrix Riccati equations [19, 20]. Matrix Riccati equations allow a "nonlinear 
superposition formula" [20], i.e. the general solution can be expressed algebraically in 
terms of a finite number of particular solutions. The superposition formula is based on 
a nonlinear action of the group SL(iV, R) with iV = 2 for the Riccati equation itself. 
A numerical method based on this group theoretical superposition formula also made 
it possible to integrate around the poles of solutions [19] and to approach the poles 
from both sides. A further relevant observation is that matrix Riccati equations can 
be discretized while preserving their superposition formulas [21, 22]. This discretization 
leads to fractional linear mappings similar in form to Equation (81) 

Table 6. Discretization errors, Example 5 



Scheme 


h = 0.1 


h = 0.01 


h = 0.001 


Sym.pres. 
Standard 


3.10 10~ 2 
5.07 10~ 3 


3.13 10" 4 
5.01 10" 5 


2.96 10~ 6 
6.70 10- 7 



6. Conclusions 

The basic motivation for this research program is that symmetries of a physical problem 
are an essential feature of the problem and should be incorporated in any mathematical 
model. In continuous descriptions, based on differential equations, this is taken for 
granted. In discrete descriptions, using difference equations, continuous symmetries 
are usually lost. It has been shown earlier [1, 2, 4, 5, 6, 7, 8, 9] that it is possible 
to construct difference schemes that possess the same symmetries as their continuous 
limits. To achieve this, it is necessary to use difference schemes (equations and meshes) 
constructed out of the invariants of the corresponding Lie groups. 

In this article, we have considered second and third order ordinary differential 
equations with three, or four- dimensional symmetry groups. Our numerical experiments 
have shown that the accuracy of the symmetry preserving schemes is much better 
(sometimes three orders of magnitude better) than that of standard schemes at no 
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significant additional cost. Example 5 has also shown that symmetry preserving schemes 
can also provide solutions when standard methods fail because of singularities. 

Imposing that symmetries be preserved in a difference scheme usually still leaves 
some freedom in the scheme. For one, or two-dimensional symmetry groups standard 
schemes are very often among the symmetry preserving ones. Starting from dimension 
three this is usually not the case. In particular all examples treated in this article are 
such that standard schemes violate the symmetries. 

We find the presented numerical experiments extremely encouraging. Future plans 
include an investigation of higher order ODEs and of systems of nonlinear ODEs from 
the point of view symmetry preserving discretizations. Also under study is the question 
of further optimizing the symmetry preserving schemes and further increasing their 
accuracy by exploiting the remaining freedom in the choice of lattices. The behaviour 
of solutions with singularities will be further studied. Finally, we are investigating 
the numerical implications of using symmetry preserving discretizations of partial 
differential equations [1, 12, 13, 14, 15]. 
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Figure 1 Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 1. 

Figure 2 Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 2. 

Figure 3 Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 3, for the case with blow-up. 

Figure 4 Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 3, for the case without blow-up. 

Figure 5 Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 4. 

Figure 6 Solution for the symmetry preserving scheme and the standard scheme, h=0.1, 
Example 5, on [0, 6]. 

Figure 7 Solution for the symmetry preserving scheme, h=0.1, 0.01, 0.001. Example 
5, on [0, 6]. 
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Figure 1. Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 1. 
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Figure 2. Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 2, h = 1. 




Figure 3. Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 3, for the case with blow-up. 




Figure 4. Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 3, for the case without blow-up. 




Figure 5. Discretization errors for the symmetry preserving scheme and the standard 
scheme, Example 4. 



Difference schemes with point symmetries and their numerical tests 



25 



500 



-500 



-1000 



-1500 



-2000 



v 



3 
X 




reference 

^v^sym.pres. 
o standard 



AAAAAAAAAAA/VWW 



Figure 6. Solution for the symmetry preserving scheme and the standard scheme, 
h=0.1, Example 5, on [0, 6]. 
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Figure 7. Solution for the symmetry preserving scheme, h=0.1, 0.01, 0.001. Example 
5, on [0, 6]. 



